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Protein molecules often self-assemble by means of non-covalent physical bonds to form extended 
filaments, such as amyloids, F-actin, intermediate filaments and many others. The kinetics of 
filament growth is limited by the disassembly rate, at which inter-protein bonds break due to the 
thermal motion. Existing models often assume that the thermal dissociation of subunits occurs 
uniformly along the filament, or even preferentially in the middle, while the well-known propensity 
of F-actin to depolymerize from one end is mediated by chemical factors (ADP complexation). Here 
we show for a very general (and generic) model, using Brownian dynamics simulations and theory, 
that the breakup location along the filament is strongly controlled by the asymmetry of the binding 
force about the minimum, as well as by the bending stiffness of the filament. We provide the basic 
connection between the features of the interaction potential between subunits and the breakup 
topology. With central-force (that is, fully flexible) bonds the breakup rate is always maximum in 
the middle of the chain, whereas for semiflexible or stiff filaments this rate is either a minimum 
in the middle or flat. The emerging framework provides a unifying understanding of biopolymer 
fragmentation and depolymerization, and recovers earlier results in its different limits. 

PACS numbers: 87.10.Mn, 82.35.Pq, 87.15.Fh 


I. INTRODUCTION 

The configuration-mediated directionality of non- 
covalent bonds between proteins explains their propen¬ 
sity to self-assemble into fibrils and filaments [TH5]. Pro¬ 
tein filaments are ubiquitous in biology, forming inside 
the cells or in the extra-cellular matrix - individually, in 
bundles, or in randomly crosslinked networks. They facil¬ 
itate the propulsion in bacteria, they control the mechan¬ 
ical strength in cytoskeleton and the bending stiffness in 
axons, they allow positional control of organelles and pro¬ 
vide the transport routes all around the cell [11 [SHE] • In a 
different situation, the self-assembly of proteins into amy¬ 
loid fibrils impairs physiological activity and is the root 
cause of a number of organic dysfunctions [51 [TTHT5] . In 
yet another context, filaments are artificially or sponta¬ 
neously assembled to achieve a specific function in the 
material, such as directed conductivity, plasmonic reso¬ 
nances, or just the mechanical strength in a fiber com¬ 
posite, with important technological applications mui]. 
Finally, a conceptually related issue emerges in the denat- 
uration of DNA [16], for which the available theoretical 
framework HZIIIH] cannot provide predictions about the 
topology of the disassembly process. The typical size of 
all these aggregates, and its time-evolution, are a non¬ 
trivial function of the rate at which bonds along the fil¬ 
ament spontaneously dissociate due to the thermal mo¬ 
tion of the assembled molecules. The dissociation rate 
and the distribution of fragments are important parame¬ 
ters which enter the master kinetic equation description 
of self-assembling filament size and populations. 

A filament growth can be summarized by the reversible 
reaction: A„ -|- Ai ^ A„_|_i, where the monomer subunit 


Ai is added to an existing filament of n-units long. For 
the forward reaction, it is commonly accepted that as¬ 
sociation proceeds by the addition of a single subunit - 
as opposed to the joining of larger segments - because 
of the greater abundance of monomers with respect to 
active fragments. In contrast, despite the importance 
of thermal breakup in many fields of colloid science and 
technology B EOj, its basic understanding is far from 
satisfactory. Several studies aimed to explain thermally- 
activated filament breakup in physical terms, came to the 
conclusion that fibrils of any respective size can aggre¬ 
gate, while the filament breakup can occur with equal 
probability anywhere along its length. In particular, 
Lee m has demonstrated that the thermal breakup oc¬ 
curs randomly along the chain, leading to daughter frag¬ 
ments of any size. In yet another classical model based 
on equilibrium detailed-balance between the various ag¬ 
gregation and breakup events, by Hill [20], the highest 
breakup probability is for two fragments of equal size, 
i.e. the breakup rate is maximum in the middle. 

Theoretical models in the past have focused on the 
simplified case of chains of harmonically bonded parti¬ 
cles (subunits), so that the binding force is linear in the 
inter-protein displacement [191120j . In this approxima¬ 
tion the normal modes of vibration of the chain are de¬ 
coupled, which makes the problem amenable to simpler 
analysis. Even in this case, previous theoretical mod¬ 
els reached contradictory conclusions, with either flat 
breakup distribution or a pronounced maximum in the 
middle. However, the physical bonds linking protein fil¬ 
ament subunits (such as hydrogen bonds and hydropho¬ 
bic attraction) are strongly anharmonic. Then the prob¬ 
lem becomes one of coupled nonlinear oscillators as in 
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the famous Fermi-Pasta-Ulam problem [^, for which 
the typical vibration modes are no longer delocalized pe¬ 
riodic waves but solitons [55]. This is also consistent 
with the finding [^|53| that in a strained Lennard-Jones 
chain, the strain is not uniformly distributed, but local¬ 
ized around the bond which is going to break first. The 
standard tools of chemical dynamics and stochastic rate 
theory |251126] , all based on the harmonic approximation 
and on normal modes, are therefore inapplicable |571|5H|. 

Here we develop a systematic microscopic understand¬ 
ing of this process based on Brownian dynamics simula¬ 
tion and theoretical arguments, focusing on the nonequi¬ 
librium breakup phenomena. Hence we study the in¬ 
trinsic breakup rates independent of any recombination 
phenomena which may occur at later stages leading even¬ 
tually to an equilibrium size. First of all, we discover 
that the topology of filament breakup critically depends 
on the bending stiffness of the chain. Secondly, a clear 
connection is found between the anharmonicity of sub¬ 
unit interaction and the fragment distribution resulting 
from thermal breakup. The anharmonic Lennard-Jones 
or Morse-like binding potential in stiff or semiflexible fila¬ 
ments inevitably leads to a very strong preference for the 
breakup to occur at chain ends, but recover the uniform, 
flat fragment distribution in the limit of harmonic (or 
any other symmetric) potential. Importantly, it is not 
the bare anharmonicity which controls this effect, but, 
more precisely, the asymmetry of the bonding potential 
about the minimum (larger force for bond compression 
than for extension), which is inherent to the most com¬ 
mon anharmonic potentials. As we will show below, it is 
precisely the asymmetry in the potential which ’’breaks 
the symmetry” between dissociation rates at the middle 
of the filament and at the ends. Those rates are equal 
only for symmetric potentials like harmonic, and they 
always differ for asymmetric potentials. 

In contrast, when the intermolecular interaction is 
purely of the central-force type, i.e. a fully flexible 
chain with no bending resistance, a bell-like distribution 
peaked in the middle is obtained in accord with the pre¬ 
diction of the Hill model. These findings can be under¬ 
stood with an argument based on counting the degrees of 
freedom per particle for the different potentials. These 
results provide a fundamental link between the features of 
intermolecular interaction and the filament breakup rate 
and topology, and can be used in the future to predict, 
control and manipulate the filament length distribution 
in a variety of self-assembly processes in biological and 
nanomaterials applications. 


II. SIMULATIONS 


To model a non-covalently bonded filament we use a 
coarse-grained model of linear chains of Brownian parti¬ 
cles (Figj^) bonded by the truncated-shifted Lennard- 


Jones (LJ) potential, 

Ulj ^ f 4g[(g/r)^^ - [cr/rf] - Uc, for r < /I'l 
ksT \ 0, for r > i?c ^ 


where r is the distance between two neighbor proteins 
i and z -I- 1, cr is the linear size of the monomer unit, 
and Uc = 4£[(cr/i?c)^^ — /Rc)^\- The parameter e = 

e/{4{{(j/R cY'^ — (cr/i?c)®] + 1} is set to maintain a con¬ 
stant well depth equal to e, independently of Re- The LJ 
potential is inherently anharmonic, except in the close 
proximity of its minimum. An alternative could be the 
Morse potential, and we have checked that the results 
do not change qualitatively with its use. Figure ex¬ 
plains what we mean by truncation: the attractive region 
stretches up to a distance Rc (indicated by arrows in the 
plot and measured in terms of LJ length scale a), while 
the depth of the potential well is kept independently fixed 
(measured by e, in units of fc^T). The shorter the attrac¬ 
tion range, the closer is the potential to its harmonic ap¬ 
proximation. For all the data we use e = 10, which well 
approximates the strength of the most common physical 
interactions such as hydrogen bonds and hydrophobic at¬ 
traction. 

We also include in our analysis the local bending en¬ 
ergy, in the form ^ K Of, where 9i is the angle between 
the directions of bonds from the particle i to the preced¬ 
ing (z —1) and the subsequent (z-l-l) subunits. Figure[^ 
illustrates the way this effect is implemented by imposing 
pairs of equal and opposite forces on the joining bonds, 
providing a net torque on the junction. It is the same al¬ 
gorithm that is used in, e.g. LAMMPS ‘angle-harmonic’ 
system [59]. The bending modulus K, in units of ksT, is 
directly related to the persistence length of the filament 
via the standard expression Ip « Ka/ksT. 

The dynamics of the chain of subunits is governed by 
the overdamped Langevin equation, 

7 ^ =-VH(r) + A(t) (2) 

where r is the vector containing the positions of all 
molecules, 7 is the friction coefficient, the total poten¬ 
tial force acting on a given particle, —W, has contribu¬ 
tion from both the LJ and the bending couples, and the 
Gaussian stochastic force defined such that (A(t)) = 0 
and {Ai{t)Aj{t')) = 2kBTj SijS{t — t') , according to the 
fluctuation-dissipation theorem. For numerical integra¬ 
tion Eq. ([^ is discretised in the form known as the 
Ermak-McCammon equation |30l|3l]: 


r{t + At) = r{t) 


VU(r) 
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At + r 



(3) 


where F is randomly extracted from a normal distribu¬ 
tion with zero average and unit standard deviation. The 
discrete time step is taken as At = 5T0“^r, where the re¬ 
duced time uint is defined as t = a'^fD, and D = k^T/j 
is the diffusion coefficient. For a typical globular protein 
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FIG. 1: (a) Scheme of the coarse-grained nanofilament as a 
sequence of subunits bonded by the truncated LJ potential, 
(b) The LJ pair interaction potential between two bonded 
subunits of size a in the chain, for several values of attrac¬ 
tive range, measured by Kc and indicated by arrows in the 
plot, (c) The contrast between a combined potential W{r) 
felt by an inner subunit in the filament, bonded on both sides, 
and the end-subunit bonded by the regular LJ potential, (d) 
Scheme of the bond-bending force which opposes changes in 
the angle between two adjacent bonds by applying couples on 
each adjacent bond. 


(e.g. Lysozyme), with diametercr ~ 5 nm and diffusion 
coefficient D ~ 10“^° m^/s [52], we obtain t ~ 0.25 /xs. 
Therefore At ~ 13 ps. Each run is initialized with the 
equilibrium interparticle distance — r^+i = 2^/®cr, as 
a straight chain (all 9i = 0), corresponding to the min¬ 
imum of all interaction potentials. A dissociation event 
is assumed to take place when one of the bonds exceeds 
the cut-off length {Rc), i.e. — r^+il > i?c, at which 
point the simulation is terminated and the location of 
the rupture recorded. The location of the rupture is 
recorded. To generate the probability distributions plot¬ 
ted in Figs. I^and]^ N independent runs are performed 
and the normalised breakup probability is calculated as 
P{s) = N{s)/N where N{s) is the total number of 
recorded breakup events for the bond s S l...n. For most 
data we have reached N > 10"^; since the runs are inde¬ 
pendent, the N(s) are binomially (Bernoulli) distributed 
and the error bars are estimated as y/P(s)[l — P(s)]/A, 
which always stayed below 10% of the value for P{s). 


III. RESULTS 

Breakup statistics along the filament 

Figure shows the main result of our Brownian dy¬ 
namics simulation: on increasing the bending stiffness 
of the filament, the highly inhomogeneous normalized 
probability P{s) changes from a bell-shaped distribution 
peaked in the middle (reminding of the Hill model), to a 
completely opposite shape, with a strong preference for 



FIG. 2: An illustration of the role of bond-bending in the 
potential. The chains of n = 20 subunits bonded with e = 10 
and Ac = 4 were initialized from a straight conformation and 
allowed to fluctuate for 5000 ts. The resulting snapshots, for 
each value of K indicated on the image, show the effect of 
differing persistence length. 



FIG. 3: The normalized probability of the first breakup as a 
function of the position along the filament for a chain n = 20 
and LJ parameters e = 10, Rc ~ 4. The effect of changing 
bending stiffness (increasing persistence length) is evident: 
for the chain with essentially no bending penalty (lowest K) 
the distribution of fragment sizes is bell-shaped with a maxi¬ 
mum in the middle of the middle of the filament. Stiff chains, 
with a strong bond-bending penalty, instead, feature a nearly 
homogeneous, flat distribution of fragment sizes - with an in¬ 
creasingly large increase of breakup rate at the ends. There is 
a broad range of semiflexible filaments that behave in exactly 
the same way: as “stiff” chains. 

single subunits to dissociate from the ends. 

The conclusion arising from this data is clear: there is 
a broad range of what one could collectively interpret as 
‘stiff’ filaments, for which the nature of bond-breaking 
statistics is exactly the same. These are with the bend¬ 
ing stiffness of AT > 1000, and their behaviour does not 
differ from the last dataset in Fig. [^(labelled ‘stiff’), cor¬ 
responding to the strictly 1-dimensional filament where 
only the motions along the chain were permitted. For 
these stiff or semiflexible filaments there is a very strong 
preference to dissociate a single subunit from the chain 
ends, which does diminish for less symmetric potentials, 
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FIG. 4: The normalized probability of the first breakup event 
as a function of the position along the filament for a chain 
n = 20 long, bonded by potentials with e = 10 and K = 1000. 
The harmonic potential with the same depth has a uniform 
(homogeneous) probability of breakup, while at increasing 
Rc (and the degree of potential asymmetry about the min¬ 
imum) the ends of the chain are increasingly more prone to 
single-particle depolymerization. The ratio of the probability 
of breaking at the end (dissociation) to the fragmentation in 
the middle Pend/Pmid = 5.46 for Rc = 4. 

as demonstrated by Fig. below. However, as the chain 
becomes increasingly flexible, the ratio of breaking rates 
at the ends and in the middle gradually reverses, and for 
a very flexible chain {K = 0.1 in the plot) the breakup 
probability resembles the prediction from the Hill model. 
One can qualitatively understand this effect: for a stiff 
filament (as shown in Figs. and[^, in order to develop 
a thermal fluctuation large enough to stretch a bond be¬ 
yond Rc, a whole sequence of bonded particles has to 
move in a correlated fashion; this leads to an effectively 
harmonic potential acting on the middle particles, and 
diminishes their breaking rate very significantly. On the 
other hand, as the particles in a flexible chain are free to 
move perpendicular to the bond axis, this coherent mo¬ 
tion is not required and the bond breaking statistics is 
dominated by the single-bond equilibrium. 

Most protein filaments are quite stiff. The F-actin has 
the quoted persistence length Ip ^ 16 /rm [331 [31], and the 
insulin amyloid filaments: Ip ^ 4/im [35] . Interestingly, 
if one measures Ip in the units of constituent protein size 
(as the parameter a in our case), these very different 
filaments all have Ip between 3000 and 6000 units. We 
therefore choose the bending stiffness modulus K = 1000 
in all subsequent analysis, which is within the class of 
‘stiff’ filaments according to the data in Fig. [^ 

This distribution of breaking points along the chain is 
equivalent to the distribution of fragment sizes resulting 
upon breakup. Figure [^ shows how this distribution de¬ 
pends on the nature of physical bond between subunits. 
As we have seen in the illustration. Fig. changing the 
cutoff distance Rc while keeping the depth of the attrac- 



FIG. 5: Relative probability of the first breakup event upon 
varying the total length N, and the parameters of bonding 
potentials e = 10, Rc = 2.5 and K = 1000. The range of en¬ 
hanced breakup probability As at each end remains constant 
for all n. 


tive potential well constant (e) effectively alters the de¬ 
gree of potential asymmetry: the larger the Rc, the more 
asymmetric the potential is. We have also independently 
tested the breaking statistics in an explicitly harmonic 
potential of the same depth and curvature at the mini¬ 
mum. In the limit of harmonic chain, we recover a com¬ 
pletely uniform (flat) distribution of fragments, with a 
very high accuracy. This is in agreement with the theory 
of Lee [ini, who assumed harmonic bonds. On the other 
hand. Fig. [^clearly demonstrates that, with increasing 
asymmetry, the breakup probability P{s) displays an in¬ 
creasingly strong preference for depolymerization from 
the ends. For the highly asymmetric (and also highly 
anharmonic) potential with Rc = 4, the breakup proba¬ 
bility of the outer bonds is over 5 times larger than the 
one of the innermost bonds. 

Another important result is shown in Fig[^ where for 
a given level of LJ potential depth and asymmetry, and 
stiff filament with K = 1000, as usual, we study the effect 
of filament length (the total number of bonded subunits, 
n). It is more difficult to normalise the breakup proba¬ 
bility P{s) this time, because for longer filaments there 
are more and more ‘plateau values’ of the constant (low) 
breakup frequency in the middle, which participate in 
the original normalisation by the total number of runs, 
N (effectively uniformly suppressing the values of P(s) 
and thus masking the characteristic ratio Pend/Pmid)- We 
therefore chose to scale all datasets by their maximum 
value of Pend, such that the different curves are compa¬ 
rable. It is clear that with the filament increasing past 
n = 20 there is no further change in the characteris¬ 
tic ratio Pend/^mid “ simply the region of ‘chain middle’ 
becomes extended. Perhaps one may regard this as an ef¬ 
fective confirmation of the Lee model [19] . since for very 
long and very stiff filaments a very large middle section 
has an effectively harmonic bonding, and therefore uni¬ 
form breakup rate. It appears, the range of enhanced 
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probability near the ends is relatively constant, As ;< 10. 
Shorter filaments have the middle region elevated simply 
because the two end-effects start overlapping. 

The finding that, for stiff filaments with asymmetric in¬ 
teraction potentials, the dissociation rate at the end can 
be substantially larger than the rate of fragmentation in 
the middle, may be important in the self-assembly kinet¬ 
ics of actin filaments mm- There, and in many other 
cases, the tendency to depolymerize at the end is ampli¬ 
fied by the presence of multiple bonds in the interior of 
the filament, due to the double-stranded helical structure 
in the case of actin. 


IV. PROBABILITY OF FIRST BREAKUP 



In addition, we studied the probability of the first 
breakup (irrespective of its position along the filament), 
upon varying Rc and the filament length n, for the case of 
a stiff filament (limit of large K) which also approximates 
the case of a ID aggregate. From the results plotted in 
log-linear fashion in Figj^it is clear that the probability 
for the chain to fracture depends exponentially on time, 
P{t) = const • with a characteristic breakup time 
increasing upon increasing the attraction range Rc 
(and the asymmetry of bonding potential with that). The 
average first-breakup time, irrespective of the location on 
the chain, is defined by = /p P(T)rdT, upon nor¬ 
malizing P{t) = A • This exponential dependence 

can be understood from the analysis of the many-particle 
Fokker-Planck equation. 


dp{r,t) 

dt 


LsPir,t) 


(4) 


with the Smoluchowski operator defined as [40l |4T] 


FIG. 6: The probability of first breakup of a filament of fixed 
n = 20, normalized such that it is equal to unity at t = 0, 
is plotted against simulation time measured in timesteps (ts). 
Different data sets represent the different attraction range 
Rc, which is our measure of potential asymmetry. The fitted 
lines are all simple exponentials, from which we extract the 
characteristic rate of the first breakup, A. 



L«(...) = DVr • [Vr(...) + f3S/ULj{r){...)] (5) 

acting on the many-particle probability density p{v,t), 
where, in supervector notation, r = {ri, ...,r„} is the set 
of interparticle coordinates. The probability as a function 
of time that all bonds remain within the cutoff at a time 
t, that is, the probability that the chain does not break 
within a time t, is given by Q{t) = p{r,t)dr. We 
shall recall that, in supervector notation, the condition 
r < Rc means that all bond vectors (relative particle 
coordinates) in the chain are within an extension smaller 
than the cutoff Rc- Furthermore, Ul.j{y) represents the 
multi-dimensional potential energy landscape given by 
the superposition of the Lennard-Jones potentials acting 
on pairs of molecules. 

The first passage/breakup time probability density 
is defined as the change of Q{t) between the time t 
and t + dt, and is thus given by P{t) = —dQ{t)/dt. 
Combining these equations, with some manipulations 
(see e.g. Ref. [42]), it is possible to show that the 
first-passage time probability density is exactly equal 
to P{t) = —DWrp{r,t)\i^ . The mean first-breakup 


FIG. 7: The mean time of the filament breakup is plot¬ 
ted for different filament lengths, indicating an almost lin¬ 
ear increase. In simple terms, taking from Fig[^ the prob¬ 
ability to break in the middle as Pmid ~ 0.035 and the 
one at the end as pend ~ 0.145, with As « 6 subunits 
from each end affected, the total rate can be estimated as 
A = (n—2As)pmid+2Aspend, whlch is the solid line in the plot 
with only a single fitting normalisation factor: 6.3 • 10“^ts“^; 
the deviations at small n are clearly due to the overlapping 
end effects (see Fig[^. 


time is then defined as the first moment of the first- 
breakup time probability density, A“^ = t ■ P{t)dt = 
t ■ [—iAVrp(r, t)|R^ ]dt, which is the same quantity as 
measured from the exponential fits in simulations. The 
exponential dependence on time can be understood from 
the analysis of the many-particle Fokker-Planck equa¬ 
tion, Eqs. Q -(§. Its general solution is p{r,t) = 
where p labels the eigenfunctions (pp 

and eigenvalues Ap of the many-body operator Lg. 
According to the ground-state dominance principle, 
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the time evolution for long filaments (n 3> 1) is dom¬ 
inated by the smallest non-zero eigenvalue Amin, such 
that, recalling the expression for P{t), the time de¬ 
pendence of the first-breakup probability is given by 
P{t) = —I?Vr/o(r, ~ g-Amint Hence the breakup 

probability is indeed exponential in time with a charac¬ 
teristic frequency-scale given by the smallest finite eigen¬ 
value Amin of the many-body operator Lg. This result ex¬ 
plains the exponential dependence on time of the breakup 
probability observed in the simulations in Fig. Also, 
combining the expressions for P{t) and for A“^, it is pos¬ 
sible to show that A « Amin, which confirms that the 
ground-state of the many-body Fokker-Planck equation 
indeed sets the time scale of breakup. 

Furthermore, the rate A grows roughly linearly with 
the chain length n, which is demonstrated in Fig. 
This particular dependence A oc n arises because the 
number of escape attempts increases with the chain size. 
One can show by means of the standard supersymmet¬ 
ric transformation of the Fokker-Planck equation into the 
Schrodinger equation [32], that A(n) is analogous to the 
quantum ground-state energy of an ensemble of (n — 1) 
bound states, and the ground state energy is extensive 
(oc n) within the quasiparticle approximation [15] . 


V. DISCUSSION 

‘Phase diagram’ of first breakup locations 



Cut-off distance Rc (a) 


FIG. 8: Contour plot showing the ratio Pend/Pmiddie as a 
function of filament stiffness K and the asymmetry parame¬ 
ter Rc- The bottom-left lagoon (dark) represents conditions 
where filaments break in the middle (bell-shaped distribution, 
according to Hill [20]), while the upper-right (light) region of 
this map represents conditions where filaments break at the 
ends (the U-shaped distribution) and negligibly in the inner 
locations. Arrows on the top signify that these geodesic lines 
extrapolate towards A —>■ oo. Arrows to the right indicate 
that there is little further change past Rc — 4 — 5. The dashed 
geodesic line marks the Pend/Pmiddie = 1 condition, separat¬ 
ing the regions of bell- and U-shape distributions. 


We find a useful representation in a map that covers all 
of the K-Rc parameter space to study how the location 
of first-breakup events along the filament changes upon 
varying both the stiffness K and the cutoff or asymme¬ 
try Rc- The results can be represented as a contour plot 
for the ratio Pend/Pmiddie as a function of K and Re- 
The contour plot is shown in Fig. [^ The bottom left 
corner, corresponding to flexible (low-A) filaments with 
short-ranged potential close to harmonic (low-i?c), repre¬ 
sents conditions where the filament breaks in the middle 
and the fragment distribution is bell-shaped, in confor¬ 
mity with Hill’s model predictions. Upon increasing both 
K and Rc at the same time, breakup in the middle be¬ 
comes less favourable and the distribution tends to flatten 
out. Eventually, for very stiff filaments and asymmetric 
potentials with large Rc the opposite limit of U-shaped 
fragment distributions with preferential breakup at the 
filament ends is recovered. This occurs in the top-left 
region of the map in Fig. [^ For symmetric binding po¬ 
tentials close to harmonic (low Re', along the K axis of 
the contour plot), the bell-shaped distribution persists 
longer upon increasing AT, eventually transforming into 
a flat distribution Pend/Pmiddie = 1 for stiff filaments. 
On the other side of the map, where Rc is increased for 
flexible chains, the bell-shaped distribution persists for 
flexible chains up to i?c oo which corresponds to the 
LJ with no cutoff. 

In general, the most dramatic change in the breakup lo¬ 


cation and fragment-distribution shape occurs along the 
path of steepest ascent, defined as the path parallel to 
the gradient of the surface. Based on our results, the 
path of steepest ascent and most dramatic evolution in 
the breakup topology is approximately identified by the 
line \og{K/kBT) = (7/5)i?c/cr. 


Bond-bending stiffness controls the filament 
breakup/recombination equilibrium 

In Figs. [^ and [^ we have shown that depending on 
the relative extent of bond-bending and central forces 
in the intermolecular interaction, the fragment size dis¬ 
tribution can change from a U-shaped distribution in 
the limit of large bond-bending rigidity, to a bell-shaped 
distribution with opposite curvature in the limit of a 
purely central-force Lennard-Jones potential. Intermedi¬ 
ate bending stiffness values yield distributions with shape 
in between the two limiting cases. 

It is first important to understand the microscopic ori¬ 
gin of this qualitative difference upon varying the bending 
stiffness in the intermolecular interaction. Since the flex¬ 
ible chain breakup statistics closely resembles the predic¬ 
tion of Hill [20] , we take a similar approach and consider 
the fragment-size dependence of the breakup rate within 
a chemical equilibrium assumption and for the special 
simplifying case of harmonic bonds. We have checked 
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that with harmonic bonds the same behaviour trend as 
in Fig. [^is reproduced, with the only difference that the 
P{s) distribution for the stiff filament is flat (as indeed 
proven by Lee [TS]) instead of U-shaped in the case of 
stiff filaments (as the last curve in Fig. shows). That 
is, the Hill-like bell-shaped P{s) is the universal result 
for fully flexible chains. 

The equilibrium constant for a dissociation reaction 
n ^ ni -f 712 of a filament n into two fragments 
rii, 772 takes the form: i^Teq = V~^Z{ni)Z{n 2 )/Z{n) = 
^nin 2 /^nin 2 ^ where Z(ni) is the partition function 
of fragment ni. K~^ „2 is the dissociation rate, while 
^ni n 2 is tiio recombination rate of these two fragments. 
The latter can be estimated from the diffusion-controlled 
collision rate of two linear chains, upon accounting for 
the diffusion coefficient of the two chains (Kirkwood- 
Riseman approximation |36] 1 and for the encounter ef¬ 
ficiency of end-to-end collisions of the two chains. In 
this way, the size-dependence was found to be ^2 oc 
(772 ln77i + 77i In772)/77i 772 77 [20j . The size-depeudeuce 
of the dissociation rate (and hence the fragment size- 
distribution) can be obtained by replacing this form for 
the association rate in the expression for the equilibrium 
constant, and upon evaluating the fragment-size depen¬ 
dence of the partition functions in the numerator of K^q. 

From classical statistical mechanics, rigid-body trans¬ 
lational degrees of freedom of the chain contribute to the 
partition function a factor ~ and rigid-body rota¬ 
tional degrees of freedom contribute an extra factor ^ , 

since the overall mass of the filament is oc 71 . Together 
these two factors give a partition function ~ 77 ®/^. The 
vibrational contributions of the monomers in the chain 
factorise in the partition function, as for a chain of har¬ 
monic oscillators, resulting in standard factors of the type 
^ (fcflT//ku)”, where to is the Einstein frequency. Clearly 
these factors do not contribute to K^q because the corre¬ 
sponding terms in the numerator and denominator can¬ 
cel. 

A full consideration of the normal modes of the linear 
chain with free ends, beyond the Einstein model, leads 
to an additional nontrivial size-dependence ~ for 

vibrations of harmonic spheres in ID, and to ^ 77 “^/^ for 
vibrations in a flexible 3D chain isaiiH]. In simple terms, 
upon increasing the chain length, more low-energy modes 
can be accommodated in the spectrum, which causes the 
partition function to decrease. The importance of this ef¬ 
fect was first recognized by J. Frenkel [35] in the context 
of nucleation phenomena. Hence with purely central- 
force interaction in 3D (flexible chain) the overall con¬ 
tribution is ^ n9/2-3/2 _ A.kin to covalent bonds in 
molecular physics, the bending stiffness introduces addi¬ 
tional degrees of freedom for rotations about the bond 
symmetry axis, which then leads to an overall depen¬ 
dence ^ 77 ®/^“^ = 77 ^/^. One should note that with 
spheres and purely central-force bonds there is no such 
axis of symmetry for the rotations, and the three transla¬ 
tional degrees of freedom per particle suffice to describe 
the vibrational behavior. Including all these considera¬ 


tions, the dissociation rate will have a dependence on the 
fragment sizes given by 

^nl.«2 ~ (771772)'^“^ (772 In 77i +77iln772)/77. (6) 

The exponent x, which collects all size-dependent con¬ 
tributions of the partition function, is different depend¬ 
ing on whether the interaction is purely central-force, or 
has a bond-bending stiffness. For central forces, x = 3, 
whereas with semiflexible or stiff chains one has x = 3/2. 
The leading contribution is then ~ ( 771772 )^, with a pro¬ 
nounced bell-shape peaked in the middle for the exclu¬ 
sively central-force flexible chain, and ^ ( 771772 )° ®, lead¬ 
ing to a much flatter distribution for a chain with bond¬ 
bending penalty. The fact that the slightly U-shaped dis¬ 
tribution observed in simulations for stiff filaments is not 
recovered by this model should be attributed to the var¬ 
ious approximations (Kirkwood-Riseman for chain diffu¬ 
sion, detailed balance, etc.) involved in the model, and 
also to the harmonic approximation of independent linear 
oscillators underlying the factorization of partition func¬ 
tions. This argument, however, explains, qualitatively, 
that a flatter distribution of fragments is to be expected 
in the presence of bond-bending, due to the additional 
rotational degrees of freedom about the stiff intermolec- 
ular bond symmetry axis, which is absent with purely 
central-force interactions. 


Possible roles of electrostatics and temperature in 
amyloid fibril breakup 

We can briefly comment on the qualitative predictions 
of this model for the distribution of breakup fragments 
in realistic amyloid fibrils. Realistic intermolecular forces 
which bind proteins in amyloid fibrils crucially depend 
on both electrostatics and temperature. We shall start 
considering the role of electrostatics first. 

Electrostatic repulsion between two bound proteins in 
a filament is ubiquitous except for solutions at very high 
ionic strength. Electrostatic repulsion acts to “lift up” 
the bonding minimum, and it may also contribute an ad¬ 
ditional small energy barrier to the total interaction U, 
with a maximum Umax co-existing or competing with the 
new lifted attractive minimum. We denote the new at¬ 
tractive minimum as < e. Due to the fact that the 

electrostatic energy decreases with r, and the maximum 
is typically at r > the lifting up of the bonding min¬ 
imum by the electrostatic repulsion is not entirely com¬ 
pensated by the energy barrier (the new maximum in U). 
Hence the total energy to be overcome for the particle to 
escape from the bonding minimum is Umax — Umin < £■ 
This consideration points towards a role of electrostat¬ 
ics which promotes breakup, or at least, restructuring 
into a different morphology where the electrostatic en¬ 
ergy density is reduced. This outcome of our analysis is 
compatible with recent experimental observations where 
an increased electrostatic repulsions (e.g. at lower ionic 


strengths) is responsible for fission or scission phenom¬ 
ena of larger compact aggregates into smaller and more 
anisotropic aggregates [nm]. 

Our simulations show a crossover from a U-shaped 
fragment distribution into a bell-shaped distribution 
upon going from high values of bond-bending stiffness 
K to lower values. In our simulations, K is fixed and 
set independently of T, the latter being kept constant 
throughout at varying K. In reality, however, K and T 
may not be decoupled for a realistic model of amyloid 
fibrils. The reason is that the inter-protein bending stiff¬ 
ness K originates, microscopically, from the strength of 
/I-sheets which bind two adjacent proteins in the fibril. 
The mechanism is known: due to the planar, sheet-like, 
nature of two hydrogen-bonded /3-sheets, there is an in¬ 
trinsic bending resistance against sliding or rolling of the 
two proteins past each other. The same mechanism pro¬ 
vides bending rigidity when two surfaces bonded by many 
anchored central-force springs are displaced tangentially 
apart. Upon increasing T, the hydrogen and hydropho¬ 
bic bonds which keep the two /3-sheets together start to 
dissociate, leading to lower bending stiffness and lower 
values of K. 

Hence, based on our simulation results, we can pre¬ 
dict that the fragment distribution function of realistic 
amyloid fibrils should evolve from a U-shaped distribu¬ 
tion at low temperature T, where the /3-sheets of two 
adjacent proteins are tightly bound, into a bell-shaped 
distribution at higher T where the /3-sheet bonding be¬ 
comes looser, which makes the bending stiffness K de¬ 
crease. This prediction seems to be confirmed by prelim¬ 
inary experiments |45] . and future work using ab-initio 
simulations should focus on identifying the relationship 
between K and T, which controls the evolution of the 
fragment distribution with T. In future research it will 
be important to combine all these effects into a general 
coarse-grained approach along the lines of [46l , to 

achieve a bottom-up description of realistic filaments and 
their size evolution. 


Anharmonicity controls depolymerization from the 
ends in stiff filaments 

When the bending rigidity of the chain is high, the 
probability of spontaneous bond breaking is flat when 
the bond potential is harmonic m - yet it adopts a 
very distinct and very strongly biased U-shape when the 
/asymmetry of the potential increases (FigQ. How can 
we quantitatively explain why the asymmetry of interac¬ 
tion potential between any two bonded subunits leads to 
higher breakup rates at the chain ends, and much smaller 
breakup rates in the middle? For a high bending modulus 
one can treat the bond at the filament end as a classical 
diatomic molecule, and a subunit in the middle of the 
chain as the inner particle in a linear triatomic molecule. 
In the latter case, the combined potential W felt by the 
particle in the middle is sketched in FigjJ. 


One would be tempted to explain the difference be¬ 
tween the higher dissociation rate at the filament end 
and the lower one in the middle by referring to the overall 
lower energy (deeper potential well) felt by the particle 
in the middle sitting in the minimum of the combined 
potential W(r). Applying a Kramers-type escape-rate 
argument would then lead to an Arrhenius dependence 
of the particle on the depth of the energy well and an 
overall large difference between the two rates. However, 
such an approach cannot explain the observation that the 
rate is the same in the middle and at the end for the case 
of harmonic potential; in that case the same argument 
about W applies hence one would expect a lower rate 
in the middle, which is not observed, in agreement with 
previous calculations m- What is different in the case 
of the harmonic potential, is the fact that the asymme¬ 
try of the bonding potential is removed for the particle 
at the end of the chain (while the subunits in the mid¬ 
dle effectively experience the harmonic potential in both 
cases). 

It is in fact this asymmetry which facilitates dissoci¬ 
ation at the termini of the chain, where less resistance 
is encountered by the particle escaping outwardly from 
the bound state. In order to verify that this is indeed the 
right physics, we also run a test simulation with a quartic 
potential U oc {r — rmin)^, which is anharmonic yet fully 
symmetric about the minimum, just like the harmonic 
potential. Also in this case we found a completely flat 
distribution of fragments, as for the harmonic potential, 
which supports the proposed claim. 

It is therefore the asymmetry, in the case of anhar¬ 
monic potentials, which plays the major role in facilitat¬ 
ing the preferential bond breakup at the chain ends. The 
explanation can be found in the different values of the 
mean thermal fluctuation from the equilibrium position 
(energy minimum) for the particle sitting in the asym¬ 
metric LJ potential at the chain end, and the particle 
moving in the more symmetric combined potential W (r) 
in the middle of the chain. An analysis of the mean ther¬ 
mal fluctuation done long ago by J. Frenkel [39], shows 
that the mean thermal fluctuation of the particle feeling 
the anharmonic/asymmetric potential at the end is typi¬ 
cally larger because of the shallower slope of the potential 
in the outward direction. For the particle in the middle, 
the situation is different because the combined potential 
W(r) does not become shallower as the particle in the 
middle moves away from one of the two neighbours, due 
to the presence of the interaction with the other neigh¬ 
bour. 


VI. CONCLUSIONS 

By means of Brownian dynamics simulations, we have 
shown that thermal breakup rates and breakup topology 
of model protein filaments (and other linear nanoparti¬ 
cle aggregates) are strongly affected by the presence of 
bond-bending stiffness in the interaction between sub- 
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units, and by the degree of asymmetry of the anharmonic 
binding potential. With stiff chains bonded by inter¬ 
particle forces with anharmonicity and asymmetry of the 
potential typical for intermolecular interaction potentials 
(van der Waals, hydrophobic attraction etc), we find a 
strongly preferential breakup at the chain ends, and an 
overall U-shaped fragment distribution. In contrast, with 
purely central-force interactions between subunits, that 
is, fully flexible chains - the fragment size distribution is 
bell-shaped, with a pronounced peak in the middle (sym¬ 
metric breakup), and the lowest breakup rate is found at 
the ends of the chain. 

While the preferential breakup at the end of stiff chains 
(filament depolymerization) can be explained in terms of 
the larger thermal fluctuations at the chain-end associ¬ 
ated with potential anharmonicity/asymmetry in a per¬ 
fectly stiff quasi-lD chain model, the dramatic change 
of breakup topology upon varying the strength of bond¬ 
bending interaction is more subtle. In this case we found 
a tentative explanation upon considering the degrees of 
freedom associated with the vibrational partition func¬ 
tion of the fragments. In general, breakup into two equal 
fragments is favoured with purely central-force bonds be¬ 
cause the product of the partition functions of two frag¬ 
ments is maximised (which is intuitive if one considers 
that the classical partition function for rigid body mo¬ 
tions increases strongly with the fragment size). The 
vibrational partition function, instead, decreases with 
fragment size because more low-energy modes can be ac¬ 
commodated in longer fragments. This effect becomes 


stronger in the case of bond-bending, where the total 
number of vibrational degrees of freedom is larger due 
to the rotation axis of the stiff bond. As a result of 
this compensation between the size dependencies of the 
vibrational and rigid-body partition functions, the size- 
dependence of fragmentation rate with bond-bending is 
much weaker compared to the central-force case. 

Hence, we found some general laws which govern the 
fragmentation behavior of model linear aggregates, as a 
function of the relative importance of central-force and 
bond-bending interactions between subunits. These find¬ 
ings are important towards achieving a bottom-up con¬ 
trol over the length and time-evolution of filament popu¬ 
lations, both in biological problems (acting, amyloid fib¬ 
rils etc.) and in nanoparticle self-assembly for photonic 
applications. 
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